int toResultIdx(int row, int col)
{
	if(row<col)
	{
		return (col-1)*col/2+row;
	}
	else
	{
		return (row-1)*row/2+col;
	}
}


__kernel void test(
	__global float *data,   //D*N
	__global float *result, //(blocksize-1)*blocksize/2, for first stage
	__local float *block1,
	__local float *block2,
	unsigned int N,
	unsigned int D)
{
	int gid = get_global_id(0);
  int wid = get_group_id(0);
	int lid = get_local_id(0);
  
  int xid = 0;
  int yid = 0;
  
  int blocksize = get_local_size(0);
  int nblocks = N / blocksize;

	float diag1;
	float diag2;
	float temp;
	float a;
	float b;

	int idx;
	
  unsigned int i;
  unsigned int j;
  
  // Determine 2D id
  for(yid = 0; yid < nblocks; yid++)
    if(wid < (yid + 1) * (yid + 2) / 2)
    {
      xid = wid - yid * (yid + 1) / 2;
      break;
    }


	diag1=0;
	temp=0;
	//row major load, and compute square terms at the same time
  for(i=0; i<D; i++)
	{
    temp=data[blocksize*yid+N*i+lid];
    block1[i*blocksize+lid]=temp;
    diag1+=temp*temp;
	}
  
  if(xid != yid) { // Varies on a per workgroup basis
    diag2=0;
    temp=0;

    for(i=0; i<D; i++)
    {
      temp=data[blocksize*xid+N*i+lid];
      block2[i*blocksize+lid]=temp;
      diag2+=temp*temp;
    }
    
    barrier(CLK_LOCAL_MEM_FENCE);

    for(i=0; i<blocksize; i++)
    {
      temp=0;
      //row-major, fast
      for(int j=0; j<D; j++)
      {
        a=block2[j*blocksize+i];
        b=block1[j*blocksize+lid];
        temp+=a*b;
      }
      
      idx=toResultIdx((yid * blocksize) + lid, (xid * blocksize) + i);
      result[idx] = diag1 - 2 * temp;
    }
    
    barrier(CLK_GLOBAL_MEM_FENCE);
    
    for(i=0; i<blocksize; i++)
    {
      idx = toResultIdx((yid * blocksize) + i, (xid * blocksize) + lid);
      result[idx] += diag2;
    }
  } else {
    barrier(CLK_LOCAL_MEM_FENCE);

    for(i=0; i<lid; i++)
    {
      temp=0;
      //row-major, fast
      for(int j=0; j<D; j++)
      {
        a=block1[j*blocksize+i];
        b=block1[j*blocksize+lid];
        temp+=a*b;
      }
      
      idx=toResultIdx((yid * blocksize) + lid, (xid * blocksize) + i);
      result[idx] = diag1 - 2 * temp;
    }
    
    barrier(CLK_GLOBAL_MEM_FENCE);
    
    for(i=lid + 1; i<blocksize; i++)
    {
      idx = toResultIdx((yid * blocksize) + i, (xid * blocksize) + lid);
      result[idx] += diag1;
    }  
  }
}